% B1_extract_wrf.m 9/27/2011 Parker MacCready
%
% this extracts fields from wrf files, and organizes them as mat files by
% variable
%
% unlike previous versions of the extraction code, this saves more of the
% domain, only clipping East of -121.5 longitude
%
% this mimics the "period" convention adopted in the mm5 processing to deal
% withtimes when the grid changed size, but it not really required here
% since the grid was always the same (so far!)

clear

% *** USER DEFINED VALUES **********************
dth = 120; % time step for forcing fields (hours) [120, 12, 3, 1]
eastlimit = -121.5; % eastward limit for extracted files
period_list = {'wrf1'};
% *** END USER DEFINED VALUES ******************

g_list = {'d2';'d3'};
indir0 = '/pmraid3/parker/tools_archive/atm/';
addpath /home/parker/mexnc /home/parker/snctools;

% make the output directory
switch dth
    case 1
        outdir0 = 'extractions_hourly/';
    case 3
        outdir0 = 'extractions_three_hour/';
    case 12
        outdir0 = 'extractions_half_day/';
    case 120
        outdir0 = 'extractions_five_day/';
    otherwise
        disp('that dth is not coded for')
        break
end
if ~(exist(outdir0,'dir')==7); mkdir(outdir0); end;

for np = 1:length(period_list)
    period = period_list{np};
    disp(' ')
    disp('**********************************')
    disp(['period = ',period]);
    
    % set the time limits of the processing
    switch period
        case 'wrf1'
             ff = '2008041512'; % first forecast
             lf = '2009123112'; % last forecast
        otherwise
            % do nothing
            break
    end
    % and create a vector of times to process for a given period
    year0 = ff(1:4); mo0 = ff(5:6); day0 = ff(7:8); hr0 = ff(9:10);
    year1 = lf(1:4); mo1 = lf(5:6); day1 = lf(7:8); hr1 = lf(9:10);
    % these are the endpoints (datenum format)
    td0 = datenum(str2num(year0),str2num(mo0),str2num(day0), ...
        str2num(hr0),0,0) + .5;
    td1 = datenum(str2num(year1),str2num(mo1),str2num(day1), ...
        str2num(hr1),0,0) + 1 - (1/24);
    % don't do last hour because the way "Z_time_to_name" is written it
    % would try to get the final hour from the start of the next period.
    tdvec = td0:dth/24:td1;
    outdir = [outdir0,period,'/'];
    if ~(exist(outdir,'dir')==7); mkdir(outdir); end;
        
    % set up the variables to process (include PP even when absent)
    varlist = {'PSFC';'RAINC';'SWDOWN';'GLW';'T2';'Q2';'U10';'V10'};
    for nvar = 1:length(varlist) % start of variable loop
        tic
        var = varlist{nvar};
        varlo = lower(varlist{nvar});
        
        for gg = 1:length(g_list) % start of resolution loop
            g_tag = g_list{gg};
            % load the lat-lon grids
            fng = [indir0,'wrf',ff(1:4),'/',lower(datestr(tdvec(1),'mmm')),'/', ...
                'wrfout_',g_tag,'.',ff,'.f12.chop.nc'];
            lon = nc_varget(fng,'XLONG');
            lat = nc_varget(fng,'XLAT');
            % clip the eastern edge of the domain
            [nr,nc] = find(lon<=eastlimit);
            ncmax = max(nc);
            lon = lon(:,1:ncmax);
            lat = lat(:,1:ncmax);
            
            % initialize result arrays
            NT = length(tdvec);
            [NR,NC] = size(lon);
            outmat = NaN * ones(NT,NR,NC);
            ts1_vec = cell(NT,1);
            ts2_vec = cell(NT,1);
            % initialize error flags
            flag_file_exists = zeros(NT,1);
            
            for tt = 1:length(tdvec) % start of time loop (hours)
                % get the file name for this time
                td = tdvec(tt);
                [ts1,ts2] = Z_time_to_name(td);
                % save the parts of the file name
                ts1_vec{tt} = ts1;
                ts2_vec{tt} = ts2;
                % make a string with the year of the forecast so that we
                % know which directory to look in
                year_fc = ts1(1:4); 
                
                % set the input directory
                indir1 = [indir0,'wrf',year_fc,'/'];
                modir = lower(datestr(td,'mmm'));
                indir2 = [indir1,modir,'/'];
                
                % form the full file name, note that it changes format at
                % bit on March 7th, 2009.
                if td < datenum(2009,3,7,12,0,0)
                    fn = [indir2,'wrfout_',g_tag,'.', ...
                        ts1,ts2,'.chop.nc'];
                else
                    fn = [indir2,'wrfout_',g_tag,'.', ...
                        ts1,ts2,'.0000.chop.nc'];
                end
                
                % now extract the data from the file, and clip the domain
                if exist(fn,'file')==2
                    flag_file_exists(tt) = 1;
                    data = nc_varget(fn,var);
                    outmat(tt,:,:) = squeeze(data(:,1:ncmax));
                end
                
            end % end of tt loop
            
            % save the output to a mat-file, note that throughout the
            % processing we use the -v7.3 flag when saving, allowing for
            % much larger file size (runs on skua, but not waddle)
            save([outdir,varlo,'_',g_tag,'.mat'],'lon','lat','tdvec', ...
                'ts1_vec','ts2_vec','outmat', ...
                'flag_file_exists','-v7.3');
            dt = toc; dts = num2str(round(dt)); nhrs = num2str(length(tdvec));
            disp(['-- ',varlo,'_',g_tag, ...
                ' (',dts,' sec for ',nhrs,' files) ------']);
            
        end % end of gg loop
        
    end % end of nvar loop
    
end % end of period loop

